Pollard $\rho$$\rho$ 算法

生日悖论

根据组合数可以得到只要一个房间里有 $23$$\color{red}23$ 个人就有 $50\mathrm{\%}$$50\%$ 的概率有两个人的生日会相同。

当 $k>56$$k>\color{red}56$，$n=365$$n=365$ 时，出现两个人同一天生日的概率将大于 $99\mathrm{\%}$$\color{red}99\%$。

考虑一个问题，设置一个数据 $n$$n$，在 $[1,1000]$$[1,1000]$ 里随机选取 $i$$i$ 个数（$i=1$$i=1$ 时就是它自己），使它们之间有两个数的差值为 $k$$k$。当 $i=1$$i=1$ 时成功的概率是 $\frac{1}{1000}$$\frac{1}{1000}$，当 $i=2$$i=2$ 时成功的概率是 $\frac{1}{500}$$\frac{1}{500}$（考虑绝对值，$k_2$$k_2$ 可以取 $k_1-k$$k_1-k$ 或 $k_1+k$$k_1+k$），随着 $i$$i$ 的增大，这个概率也会增大，最后趋向于 $100\mathrm{\%}$$100\%$。

随机数

我们设 $f(x)=(x^2+c)modn$$f(x)=(x^2+c)\bmod n$，其中 $n$$n$ 是待分解的整数；

再设数列 $x$$x$ 的递推式为 $x_i=f(x_{i-1})$$x_i=f(x_{i-1})$。初值为 $x_0=1$$x_0=1$。

这样可以获得足够的随机性，并且根据生日悖论，只要生成 $\sqrt{n}$$\sqrt{n}$ 个数就可以找到环。

设 $m$$m$ 为 $n$$n$ 的最小非平凡因子，则一定有 $1<m<n$$1<m<n$。

于是我们再设数列 $y$$y$，满足 $\mathrm{\forall }i\in [1,\mathrm{\infty }],y_i=x_{i}modm$$\forall i\in\lbrack 1,\infty\rbrack,y_i=x_i\bmod m$。如果 $\mathrm{\exists }i,j$$\exists i,j$ 使得 $i\ne j,x_i\ne x_j,y_i=y_j$$i\neq j,x_i\neq x_j,y_i=y_j$，我们就找到了 $n$$n$ 的最小非平凡因子 $m$$m$。

再应用一次生日悖论可以得到时间复杂度为 $O(\sqrt{m})\le \mathrm{\Theta }(n^{\frac{1}{4}})$$O(\sqrt{m})\leq\Theta(n^\frac{1}{4})$。

Floyd 判环算法

假设两个人同时同向同地开始赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会追上 B，且 B 第一次被追上时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 $1$$1$ 倍。

于是我们可以写出这样的代码：

inline __int128 abs(__int128 x) { return x < 0 ? -x : x; }
#ifdef __USE_MY_RAND__
static unsigned long long seed;
unsigned long long u64_hasher(unsigned long long x) {
    return x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x * 0x1952BB648B74B19Eull;
}
unsigned long long rng() { return seed += u64_hasher(seed); }
#else
static std::mt19937_64 rng;
#endif
inline uint64_t rho(uint64_t value) {
    if (value % 2 == 0) return 2;
    if (value % 3 == 0) return 3;
    if (value % 5 == 0) return 5;
    if (value % 7 == 0) return 7;
    mint::set_mod(value);
    mint x, y;
    mint c;
    while (true) {
        c = rng();
        y = c;
        x = c * c + c;
        while (x != y) {
            mint z((uint64_t)abs((__int128)y.val() - (__int128)x.val()));
            if (uint64_t g = std::__gcd(z.val(), value); g > 1 && g < value) return g;
            y = y * y + c;
            x = x * x + c;
            x = x * x + c;
        }
    }
}

倍增优化

上面的做法要调用很多次 $gcd$$\gcd$，导致算法乘上了一个 $\log$$\log$。

其实这个 $\log$$\log$ 是可以避免的（和常数抵消）。

修改上述结论：假设两个人同时同向同地开始赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会追上 B，且此时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 $k(k\in {\mathbb{Z}}^{+})$$k(k\in\mathbb{Z}^+)$ 倍。（无论之前 B 有没有被 A 追上过）

我们每过一段时间将这些差值进行 $gcd$$\gcd$ 运算。设 $z=\prod |x_i-y_i|modn$$z=\prod|x_i-y_i|\bmod n$，由欧几里得算法得到 $gcd(abmodn,n)=gcd(ab,n)$$\gcd(ab\bmod n,n)=\gcd(ab,n)$，所以 $gcd(z,n)=gcd(\prod |x_i-y_i|,n)$$\gcd(z,n)=\gcd(\prod|x_i-y_i|,n)$。

如果某一时刻得到 $z=0$$z=0$ 那么表示分解失败，退出并返回 $n$$n$ 本身。每隔 $k$$k$ 个数，计算是否满足 $1<gcd(z,n)<n$$1<\gcd(z,n)<n$。此处取 $k=127$$k=127$，可以根据实际情况进行调节。

更正：由拉梅定理得到 $gcd$$\gcd$ 的辗转相除次数 $\le 5{\log }_{10}\text{min}(n,m)=5{\log }_{10}2\times \log \text{min}(n,m)$$\leq 5\log_{10}\text{min}(n,m)=5\log_{10}2\times\log\text{min}(n,m)$，所以对于 $64$$64$ 位整型可取 $k=100$$k=100$。

代码

用几个小素数判一下效率会有显著提升（？）。

附带 abs(__int128)、readu128()、putu128()。这份代码中采用倍增优化，$k=32768$$k=32768$。

这里有一个小技巧：只改变 $x$$x$ 在环上的位置，不改变 $y$$y$ 在环上的位置。可以证明这和一个前进两格、另一个前进一格是等价的。

inline __int128 abs(__int128 x) { return x < 0 ? -x : x; }
#ifdef __USE_MY_RAND__
static unsigned long long seed;
unsigned long long u64_hasher(unsigned long long x) {
    return x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x * 0x1952BB648B74B19Eull;
}
unsigned long long rng() { return seed += u64_hasher(seed); }
#else
static std::mt19937_64 rng;
#endif
inline __uint128_t rho(__uint128_t value) {
    if (value % 2 == 0) return 2;
    if (value % 3 == 0) return 3;
    if (value % 5 == 0) return 5;
    if (value % 7 == 0) return 7;
    mint::set_mod(value);
    mint x, y, z, c;
    uint64_t i, j;
    __uint128_t g;
    while (true) {
        y = x = rng();
        z = 1;
        c = rng();
        i = 0;
        j = 1;
        while (++i) {
            x = x * x + c;
            z *= mint((__uint128_t)abs((__int128)y.val() - (__int128)x.val()));
            if (x == y || !z.val()) break;
            if (!(i & 32767) || i == j) {
                g = std::__gcd(z.val(), value);
                if (g > 1ull) return g;
                if (i == j) {
                    y = x;
                    j <<= 1;
                }
            }
        }
    }
}
inline auto crack(__uint128_t value) {
    struct _RetType { __uint128_t base; uint64_t power; };
    std::vector<_RetType> ret;
    if (prime(value)) {
        ret.push_back({value, 1ul});
        return ret;
    }
    _RetType node{2, 0};
    while (value % 2 == 0) ++(node.power), value >>= 1;
    if (node.power) ret.push_back(node);
    auto dfs = [&](auto self, __uint128_t v) {
        if (v <= 1) return;
        if (!prime(v)) {
            __uint128_t tmp = rho(v);
            self(self, tmp);
            self(self, v / tmp);
        }
        else {
            auto pos = std::find_if(ret.begin(), ret.end(), [=](const auto& x) { return x.base == v; });
            if (pos == ret.end()) ret.push_back({v, 1ul});
            else ++(pos->power);
        }
    };
    dfs(dfs, value);
    std::sort(ret.begin(), ret.end(), [=](const auto& x, const auto& y) { return x.base < y.base; });
    return ret;
}
inline __uint128_t readu128() {
    __uint128_t ret = 0;
    char c = getchar();
    while (!isdigit(c)) c = getchar();
    while (isdigit(c)) ret = ret * 10 + (c & 15), c = getchar();
    return ret;
}
inline void putu128(__uint128_t x) {
    static char buf[45];
    static int curpos;
    curpos = 0;
    do {
        buf[curpos++] = x % 10;
        x /= 10;
    } while (x);
    while (curpos--) putchar(buf[curpos] | 48);
}